#ifndef OUTPUT_HPP
#define OUTPUT_HPP

#include "Multigrid_One_Dimension.hpp"
#include "Multigrid_Two_Dimension.hpp"
#include "Function.hpp"
#include "json.hpp"
#include "Input.hpp"
#include <iostream>
#include <fstream>

using json = nlohmann::json;
double output_1d(const Eigen::VectorXd &solution)
{
    // Read input file
    std::ifstream input_file("../Input.json");
    if (!input_file.is_open())
    {
        std::cerr << "Error: Could not open input file!" << std::endl;
        std::abort();
    }

    json input_data;
    input_file >> input_data;
    std::string source_term_str = input_data["source_term"].get<std::string>();
    std::shared_ptr<Function_One_Dimension> function;
    if (source_term_str == "PI^2 * sin(PI x)")
    {
        function = std::make_shared<function1>();
    }
    else if (source_term_str == "e^(x^2)")
    {
        function = std::make_shared<function2>();
    }
    else if (source_term_str == "-(-sin(x) + cos^2(x) + 1) * exp(x + sin(x))")
    {
        function = std::make_shared<function3>();
    }
    else if (source_term_str == "2 * PI^2 * sin(PI * x) * sin(PI * x)") {
        function = std::make_shared<function4>();
    }
    else if (source_term_str ==  "-(-cos(x) - sin(x) + sin^2(x) + cos^2(x)) * exp(sin(x) + cos(x))") {
        function = std::make_shared<function5>();
    }
    // 输出解
    std::cout << "Solution:" << std::endl;
    int m = solution.size() - 2;
    double h = 1.0 / (m + 1);
    // for (int i = 0; i < solution.size(); ++i)
    // {
    //     double x = (i)*h;
    //     std::cout << "x = " << x << ", u(x) = " << solution[i]
    //               << ", exact = " << function->operator()(x) << std::endl;
    // }

    // 计算误差
    double error_norm = 0.0;
    for (int i = 0; i < solution.size(); ++i)
    {
        double x = (i)*h;
        error_norm += std::pow(solution[i] - function->operator()(x), 2);
    }
    error_norm = std::sqrt(error_norm);

    int n = solution.size();
    char filename[100];
    sprintf(filename, "test_parm%d.txt", n-1); // 修正文件名生成
    FILE *fp = fopen(filename, "w+");
    if (!fp)
    {
        std::cerr << "Failed to open file: " << filename << std::endl;
    }
    // 输出数值解到文件
    for (int i = 0; i < n; i++)
    {

        fprintf(fp, "%lf %lf\n", i * h, solution[i]);
    }
    fclose(fp);
    return error_norm;
}

double output_2d(const Eigen::VectorXd &solution)
{
    // Read input file
    std::ifstream input_file("../Input.json");
    if (!input_file.is_open())
    {
        std::cerr << "Error: Could not open input file!" << std::endl;
        std::abort();
    }

    json input_data;
    input_file >> input_data;
    std::string source_term_str = input_data["source_term"].get<std::string>();
    std::shared_ptr<Function_Two_Dimension> function;
    if (source_term_str == "-(-sin(x) + cos^2(x) + 1) * exp(y + sin(x))")
    {
        function = std::make_shared<function1_2d>();
    }
    else if (source_term_str == "2 * PI^2 * sin(PI * x) * sin(PI * y)")
    {
        function = std::make_shared<function2_2d>();
    }
    else if (source_term_str == "-(-cos(y) - sin(x) + sin^2(y) + cos^2(x)) * exp(sin(x) + cos(y))")
    {
        function = std::make_shared<function3_2d>();
    }
    else
    {
        std::cerr << "Error: Invalid source term!" << std::endl;
        std::abort();
    }
    // 输出解
    std::cout << "Solution:" << std::endl;
    int m = std::sqrt(solution.size()) - 2;
    double h = 1.0 / (m + 1);
    // for (int i = 0; i < m + 2; ++i)
    // {
    //     for (int j = 0; j < m + 2; ++j)
    //     {
    //         double x = (i)*h;
    //         double y = (j)*h;
    //         int t = i + j * (m + 2);
    //         std::cout << "x = " << x << ",  y = " << y << ", u(x) = " << solution[t]
    //                   << ", exact = " << function->operator()(x, y) << std::endl;
    //     }
    // }

    // 计算误差
    double error_norm = 0.0;
    for (int j = 0; j < m + 2; ++j)
    {
        for (int i = 0; i < m + 2; ++i)
        {
            double x = (i)*h;
            double y = (j)*h;
            int t = i + j * (m + 2);
            error_norm += std::pow(solution[t] - function->operator()(x, y), 2);
        }
    }
    error_norm = std::sqrt(error_norm);

    int n = std::sqrt(solution.size()) - 1;
    char filename[100];
    sprintf(filename, "test_parm%d.txt", n); // 修正文件名生成
    FILE *fp = fopen(filename, "w+");
    if (!fp)
    {
        std::cerr << "Failed to open file: " << filename << std::endl;
    }
    // 输出数值解到文件
    for (int i = 0; i < m+2; i++)
    {
        for (int j = 0; j < m+2; j++)
        {
            int t = i + j * (m+2);
            double x = (i)*h;
            double y = (j)*h;
            fprintf(fp, "%lf %lf %lf\n", x, y, solution[t]);
        }
    }
    fclose(fp);
    return error_norm;
}

#endif // OUTPUT_HPP